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Abstract 

We describe a method to model nonlinear dynamical systems using 
periodic solutions of delay-differential equations. We show that any finite- 
time trajectory of a nonlinear dynamical system can be loaded approxi- 
mately into the initial condition of a linear delay-differential system. It 
is further shown that the initial condition can be extended to a periodic 
solution of the delay-differential system if an appropriate choice of its 
parameters is made. As a result, any finite set of trajectories of a non- 
linear dynamical system can be modeled with arbitrarily small error via 
a set of periodic solutions of a linear delay-differential equation. These 
results can be extended to some non-linear delay differential systems. One 
application of the method is for modeling memory and perception. 

1 Introduction 

How the information about the outside world is stored in the brain is still in 
many ways an open question. Clearly, representation of the totality of the in- 
formation involves some sort of modeling of the environment by the brain under 
the specific conditions of its operation, e.g., the limited number of the neurons 
and synapses, the inherent delays in signal propagation along the neural fibers, 
etc. Assuming that the environment is described by some nonlinear determin- 
istic dynamical system and that the brain is described by another dynamical 
system, this question can be rephrased as the question of how one dynamical 
system can model another. Since the environment is typically much more com- 
plicated than the brain, the modelling task should be impossible simply on the 
count of the difference in the number of the degrees of freedom for the two 
dynamical systems. A possible solution to this puzzle is to assume that the 
environment is actually a collection of weakly interacting subsystems, each of 
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which has a lower complexity than the brain. A good strategy to learn com- 
plicated environment would then be to model the strongly coupled degrees of 
freedom group by group. Decomposition of the environment into a collection 
of components that interact weakly presumably should involve some version of 
nonlinear factor analysis extended to deterministic systems. 

Even if the complexity of the environment is less than that of the system 
that models it a fundamental question arises about in what sense and how 
one dynamical system can model another. In addition one can ask whether 
there exist a class of dynamical systems that in some sense are universal in 
their modeling properties in the sense that a large class of nonlinear dynamical 
systems that represent possible environments can be modeled using essentially 
the same system. Of course, in order for this to make sense in explaining the 
properties of the brain, the modeling has to be physically realizable. 

Leaving investigation of the deterministic nonlinear factor analysis to an- 
other publication |l| , in this paper we present arguments that when the dimen- 
sion of the environment is less or equal to that of the modeling system such 
universal modeling systems exist and that it is by the essential use of the delays 
that the universality can be realized. We prove that under certain conditions 
any nonlinear dynamical system can be modeled by a linear delay-differential 
dynamical system with constant coefficients, in the sense that any finite set 
of trajectories of the environment can be loaded into periodic solutions of the 
delay-differential equation. The loading turns out to be equivalent to solving 
an essentially linear problem and is achieved by adjusting the coefficients of the 
delay-differential equation in a way that is plausible from the biological point 
of view. With the biological applications in mind we also extend some of the 
results to a class of nonlinear systems with delays. 

Modeling the environment is usually referred to as construction memory 
models of the environment. Since the early eighties much work has been done 
on memory models under the assumption that time-independent memories can 
be represented as fixed points of systems of nonlinear differential equations 
(ODE) that possess a Liapunov function and with the evolution given by 

±{t) = F(x), x(t) e R N (1) 

for some function F (x) . Once a Liapunov function is given, the fixed points 
of the evolution can be identified with the iV-dimensional stored patterns. Do- 
mains of attraction of the fixed points then can be viewed as all the patterns 
that will be " associatively" recalled dynamically as the stored pattern repre- 
sented by the fixed points. Especially useful within this approach proved to be 
the Hopfield model || described as a dynamical system with evolution 

x(t) = -E M -x(t) + A-a(x(t))+y, (2) 

where = diag{fj,i}, fa > is an N x N diagonal matrix, A - is an N X N 
matrix, called the matrix of weights, a : Xi — > <Ji(xi) is a component-wise nonlin- 
ear transformation and the constant vector y £ R N describes time-independent 
external inputs to the system. If o'(z) > and matrix A is symmetric, or if 



2 



the symmetric part of A is non-positive definite ||, then the system (^) has a 
Liapunov function V (x) defined on the configuration space such that on the 
trajectories V(x(t)) < 0, leading to the convergence of all trajectories to a set of 
fixed points. An extensive theory exists about the properties of memory storage, 
especially for random patterns [H . 

A certain amount work has been done for the systems with delays. Existence 
of a Liapunov function can be proven for the Hopfield model with delays if the 
delays are not too large [9. However, the fixed-points paradigm is difficult to 
apply to the time dependent patterns. Time-variable patterns, which in our 
approach are represented by the trajectories of nonlinear dynamical systems, 
are ubiquitous in the environment. Therefore, understanding the underlying 
principles of their storage is important when constructing plausible models for 
perception and memory. Some interesting effects of the presence of delays on the 
structure of the attractors leading to multistability in networks of two Hodgkin- 
Huxley type spiking neurons were discussed heuristically and numerically in . 
From our point of view the model that is considered there is that of a system 
of two linear delay differential equations with time variable coefficients. In this 
paper we do not restrict ourselves to the models of spiking neurons but consider 
arbitrarily large systems of linear and nonlinear delay-differential equations with 
slowly varying coefficients. Systems with time variable coefficients will be con- 
sidered in detail elsewhere. Nevertheless, our results are in general agreement 
with the numerical simulations in j(| and provide a new perspective on their 
interpretation. 

Here we present an approach to the time-dependent pattern storage that 
describes sets of interacting " neurons" as systems of delay-differential equations 
(DDE) with multiple delays. Since delays in nerve impulse propagations are 
common, the approach is more plausible biologically than, for example, the 
Hopfield model where instantaneous communication among the neurons is as- 
sumed. Mathematically the delays appear after integrating out some internal 
degrees of freedom in systems of Hodgkin-Huxley equations (tJ . Although more 
complicated to analyze, the DDE have some properties that ordinary differential 
equations do not have. For example, for a solution of a DDE to be unique, one 
has to specify an initial condition which is a function defined on a finite time 
range This is a key feature of DDE that we exploit. 

Our starting point is to include explicitly the multiple delays ti,...,tx in 
the dynamics of memory model so that instead of (|l|) we obtain a system of 
nonlinear DDE 



x(t) = F [x (t) ,x (i - n) , ...,x {t ~ TL )) . (3) 

One example of such a system is a generalization of the Hopfield model to 
include delays. This was considered in Q 

L 

x{t)=^A k -a(x(t-T k ))+y{t), (4) 

k=0 
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where as before the Ah and are N x N weight matrices with Aq = —E^, and 
the delays are defined so that r = 0, Tfc < t^+i- In contrast to @ we assume 
the external inputs y(t) as time- variable. In models describing networks of oscil- 
lating neurons one can also consider Tk as an N x N matrix with each element 
{Tk)ij describing the delay in signal propagation from neuron (J) to neuron (i) 
via "fc"th pathway. However, after appropriate redefinition of matrices A^ and 
nonlinearities a (z) one can recover (Q) at the expense of increasing L to the 
number of different values of all matrix elements of all (L + 1) matrices. We 
assume that y[t) is a trajectory of some nonlinear dynamical system with local 
evolution given by y(t) — G(y). The detailed nature of G(y) is not important 
for the discussion. When <Ji(z) — z we obtain a system of linear DDE. Writing 
it in components we obtain 

L 

^ = J2 ( Ak ^j ■ x j (* - Tfc ) + y> (*)> * = x > N - ( 5 ) 

The main point of this paper is to prove the existence of the periodic solutions 
of systems of homogeneous linear DDE that are periodic extensions of their 
initial conditions. That is, we show that, given a set of fixed delays Tk, the 
weight matrices A^ can be chosen in such a way that the initial condition when 
extended periodically to t > tl is the solution of ([5]) . With appropriate choice of 
parameters the choice of A^ is unique. We also show that if the inhomogeneous 
part of a DDE is given by a trajectory of a nonlinear dynamical system, the 
trajectory can be loaded approximately, but with arbitrarily small error into the 
initial condition for the DDE and hence stored as its periodic solution. With 
some modifications these results can be carried over to a class of nonlinear 
DDE. Using the method we describe, with a sufficiently large N, any finite 
set of trajectories of a nonlinear dynamical system, for example any finite set 
of periodic orbits, can be encoded as periodic solutions of a linear DDE. Since 
many dynamical systems, are uniquely characterized by the sets of their periodic 
orbits, we can speak of linear modelling of nonlinear systems. 

The paper is structured as follows. In the next section we give a brief 
overview of the differences between the linear ODE and the linear DDE from 
mathematical point of view. Section 3 contains discussion of loading the tra- 
jectories into the initial conditions. In Section 4 we construct solutions for the 
periodic extension problem. In Section 5 we discuss the adaptation dynamics 
for the weight matrices Ak ■ Section 6 is a summary. 

2 Delay-differential Equations 

To specify a solution of a differential equation uniquely one needs to pose the 
initial value problem. There is a big difference in how the initial value problem 
is posed for ODE and for DDE. For the ODE the initial value problem is given 

by 

x (4>, y) (i) is a solution of (Q) for t > 0; 
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x(<j>,y) (0) = (f) £ R 



For the DDE of retarded type (||) with the maximum delay tl the initial 
value problem is formulated as 

x (4>, y) (t) is a solution of (^) for t > tl] 
x (0, y) (0) = <j> (t) € R N for t G [0, r L ] . 

Hence to specify a unique solution of a DDE one needs an initial value from 
a space of functions defined on the interval [0,tx]. Exploiting the difference is 
one key feature of our approach. 

Let us consider in more detail the difference in the linear case. Similar to 
homogeneous linear ODE we can search for the exponential solutions of ho- 
mogeneous DDE. As for ODE these correspond to zeroes of the characteristic 
equation that is obtained using the Laplace transform. For linear ODE the 
characteristic equation is polynomial and has exactly N complex roots. The 
situation is different for linear DDE. The characteristic equation becomes tran- 
scendental and is given by 

det {H (A)) = 0, (6) 

L 

H(X) = A-/-^A fc exp(-Ar fc ). 

k=0 

There are generically infinite number of complex roots A p of Jkj) . Let N(a,b), 
a,b G R, a < b < oo, be the number of the complex roots of (g) in the vertical 
strip of the complex plain with real part of each point contained in [a, b] . Then, 
some of useful properties of the roots are 

(i) Re(Ap) < c < oo if all Tfc > 0, for all k,i,j; 

(ii) Re(Ap) — > — oo p — > oo; 
(hi) N(a, b) < oo. 

Solutions of linear inhomogeneous DDE can be conveniently written down in 
an integral form using the notion of the fundamental solution. The fundamental 
solution is defined as the unique matrix solution of (JsJ) with y (t) = and with 
the following initial condition 

X (t) = t g [0, tl] ; 

X (t) = I t = tl, where / is the unit matrix. 

Using inverse Laplace transform, the fundamental solution can be expressed 
in terms of the characteristic (matrix) function H (A) 

X(t) = ± [ d\ e xt H~ x (A) , (7) 

2m J(c) 
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where (c) is the contour (c — ioo, c + zoo) in the complex plane such that Re (A p ) < 
c and where the Laplace transform and its inverse are defined by 



x (A) = / dX exp (-A t) x (t) , (8) 
Jo 

x(t) = ^— I dX exp (-A t) x (A) . (9) 
2tti J (c) 



With these definitions the unique solution of the initial value problem for (||) 
can be written as [|| 

x{<t>,y) (t) =X(t-T L )<t>(r L )+ I dTX(t-r)y(r) (10) 



TL 



YsT drX(t-T-T k ) A k <f>(T). 

I 1 J TI. —Ti- 



k=1 J T L -T k 

If we put L = then the third term in (|o[) vanishes and we recover the solution 
of a linear inhomogeneous ODE in terms of its fundamental solution. Note that 
in (|lo|) the inhomogeneous part y (t) and the initial condition <fi (t) play a similar 
role. We shall exploit the similarity below when we generate initial conditions 
from the inhomogeneous part of the equation. 

3 Loading Trajectories into the Initial Condi- 
tions 

In this section we describe how to generate the initial condition for a DDE from 
its inhomogeneous part. The main idea is to consider the inhomogeneous part 
to be nonzero only on [0, tj\ and generate the initial condition recursively by 
adding more and more terms to the equation successfully. 

Let us first consider the linear DDE. We divide the interval [0, tl] into ad- 
joining intervals [t^ T fc+1 ] , k — 0, ...,L — 1 and for time t < r m truncate the 
full DDE to 

m 

x{t)=^A k -x{t-T k ). (11) 

fc=0 

This can be done under the assumption that the signals with larger delays did 
not arrive yet at the "neuron" to influence its state. 

When m = ( |TT| ) reduces to a linear ODE which can be integrated for t € 
(0, Ti] using the initial condition <^>o at t = and the values of the inhomogeneous 
part y(t). The result of the integration can be considered as the initial condition 
the truncated DDE defined on t G (ti,T2], which enables us to integrate the 
DDE with one delay Tion the segment [ti.t 2 ] using y(t) defined on [ti^] only. 
Proceeding in this fashion step by step we can construct the initial condition 
on the entire segment [0,tl]. Taking the inhomogeneous part y(t) to be zero 
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outside of [0, tl] we can interpret the result of this procedure as a linear mapping 
of trajectories into the initial conditions, II ■ y(t) — » $ (t) , t e [0, tl]. 

To describe the iteration procedure in more detail we give the initial iteration 
step followed by the "m"th iteration. For the initial step we need to generate 
the initial condition for 

x (t) = Aq ■ x (t) +A 1 -x(t-n)+ y(t), (12) 
y(t)=0,tt[0,n}. 

These are given by the solution of the initial value problem 

x(t) = A -x(t)+y(t), (13) 

x(0) = ho- 
using the fundamental solution for this equation the solution for t > and the 
initial condition for (|l^) can be written as 

$o (t) = X a (t) ■<f> + [ ds X (t - s) y(s) , (14) 



(i 



X (t) = exp (A t) 



With this initial condition the solution for ( |i~2[ ) can be written in terms of its 
fundamental solution as 

*i (t) =X 1 (t- n) $ (n) - ^ dsX x {t-s- ti) Ax $ (*) (15) 

Jo 

+ { ds Xx (t-s) y (s) . 



$i (t) can be considered as the initial condition for the 2nd step of iteration, 
the result of which is the initial condition defined on the segment [ti^] • The 
general form of iteration is easily deduced. Namely, 



$ k (t)=X k (t-T k ) $ k -x(T k )- 

+ f ds X k (t-s) y (s) . 



f dsX k {t~s~ T k ) A k $ fc _x (s) (16) 



with the final step resulting in the initial condition 

$l (t) = {$ fe (t) , t e [T k ,n+x] ,k = o, L — i}. (17) 

Next note that to enable its periodic extension the initial condition has to be 
periodic, since it has to be defined on a closed segment. Hence we obtain a 
constraint of the form 

$L (r L ) = $L (0) = 0o- (18) 
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This condition can be used to eliminate an unknown parameter 4>q and make 
trajectory loading unique. For example, when there is only one delay we obtain 

0o = (exp(-A r L ) - 1)~ / dsexp(-A s) y(s). (19) 

Jo 

When multiple delays are present 4>q can be determined similarly. The null 
space of the mapping II is also of interest. It describes all the trajectories that 
cannot be loaded into the periodic extensions of the initial conditions. For the 
example with only one delay this space consists of all functions with the Laplace 
transform satisfying 

y (A) + O = (20) 

For arbitrary set of delays the null space still forms a one-parameter family 
and, therefore, its existence excludes only a very small set of all trajectories 
from loading. 

Obviously, the detailed loading procedure described above works only for 
the linear DDE. However the same principle can be applied to the nonlinear 
DDE as well. If the delays are known, then one can proceed with the same 
iterative scheme by numerical integration. We speculate that some sort of the 
analog integration version of the procedure above might be used in the brain to 
implement the loading. With this in mind we can interpret the time tl as the 
" attention span" of the population of the neurons that encode a particular time- 
dependent pattern. For neurons with a large number of synapses and, hence a 
large number of random delays, a continuous version of the iteration in ( |l6| ) 
can be easily written. 

4 Periodic Extensions of the Initial Conditions 

Having shown how trajectories can be loaded into the initial conditions of the 
linear DDE, we now construct solutions of the DDE that are the periodic ex- 
tension of their initial conditions. First we consider the linear case and then 
discuss the modifications added by the presence of nonlinearities. 

4.1 Linear DDEs 

Consider a Fourier series representation of a space of periodic and continuous 
initial conditions on [0, t£\ 

Q 

4 (*) = J2 ex P ' 6 «' 1 e [°> T ^ ' ( 21 ) 

n=-Q 

2ir 

Pn = — n, 

where b n are fixed iV-dimensional complex-valued vectors such that under com- 
plex conjugation b n = b- n (to ensure that 4>{t) is real- valued) . If the initial 
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condition does not belong to this space we shall consider ( f2l| ) as an approxima- 
tion of the initial condition with the error of approximation determined by the 
truncation of the infinite Fourier series that represents it at the "Q"th term of 
the expansion. 

Let us extend the domain of definition of the initial condition <f> (t) to all t > 
by treating <f> (t) as a continuous periodic function for t > with tfi (t + tl) = 
<fi (t) . Assume now that all delays are fixed but we are free to vary the matrices 
Ak. The problem of finding periodic extensions is then to find such Ak that for 
t > the function <j) (t) is a solution of the homogeneous equation 

L 

j>{t) = J2 A k-<P(t-T k ). (22) 

k=0 



After the substitution of (|21|) into (J22j) we obtain a set of linear equations on 
matrices Ak 

ip n I - ^2 A k ■ cxp (-ip n Tk) •b n = 0, n = 0, Q. (23) 
V fc=0 / 

Note that the equations for n = — Q, — 1 are redundant, since Ak are real 
and hence the additional equations can be obtained from ( |2^ ) by complex con- 
jugation. Altogether, because the n = equation is real, ( p3| ) is equivalent 
to (2Q + 1) N real equations on (L + 1) iV 2 elements of Ak, provided that all 
elements of Ak are generically non-zero. In the case a delay Tk is given by a 
random value of delays in propagation from "i"th to "j" neurons, with proba- 
bility one only one of the elements of Ak is non-zero (see the discussion in the 
Introduction). However, in such a case, if L' + 1 is the number of N x N delay 
matrices, then L + 1 = (£' + 1) N 2 and, therefore, the count of unknowns is 
the same. We conclude that ( |23{ ) can have a unique solution for Ak only when 
(2Q + 1) = (L + 1) N. Since Q determines the number of terms in the Fourier 
expansion and, therefore, the accuracy of representation of an arbitrary initial 
condition, this relation means that for a given number L of delays and N the 
number of neurons involved in modeling one can achieve only limited accuracy 
of representation. Namely, the number of expansion terms Q is bounded so that 
(2Q + 1) < (L + IV N. 

The system (|23|) can be solved by noticing that it implies that vectors b n 
are eigenvectors of linear combinations of Ak with eigenvalues ip n . Unlike the 
typical eigenvalue problem where Ak are known and one needs to find p n and 
b n , here the situation is reversed: one needs to find Ak assuming that p n and 
b n are known. If an eigenvalue A and an eigenvector b are given then all the 
matrices D for which they solve the eigenvalue problem can be written as 

d = xi + b(i- \\b\r 2 b® b 



where B is arbitrary N x N matrix and ||6|| 2 = Xa=i l^*| 2 an d b®b is the matrix 
that is the outer product of b with its complex conjugate. For given eigenvalue 
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and eigenvector matrix D has N 2 — N free parameters. Therefore, an equivalent 
way to write (123) is 



Y^Rnk -A k =i Pn I+B n (j- ||6„|r 2 6 n ®5„) , n = 0,...,Q, (24) 



k=0 



i? nfe = exp {-ip„T k ) 



(25) 



where -B n is an N x iV complex matrix, which effectively has iV 2 — N free 
parameters. This linear system of equations always has solutions for appropriate 
choice of parameters Q, L, N. For some choices the solution is unique. This 
concludes the proof of existence of the periodic extensions. 

Let us consider ( 23} 24j) for various choices of parameters. When N = 1 
the second term in the RHS of (24) vanishes and the dependence on b n drops 
out. Consequently, the system (23) becomes a system of (2Q + 1) real linear 
equations on L + 1 real parameters A k 



^ SnkAk 



,Q, 



k=0 
L 



J2 C nkA k = 0, n = 0,l,...,Q, 



fc=0 



Snk = Sin (p n Tfe) . 
C nk = COS (p n T k ) 



(26) 

(27) 

(28) 
(29) 



When L = 2Q, provided detT ^ 0, the equations can be solved uniquely by 
inversion of the matrix T 



T = 



( sin (piT ) 

sin (pqt ) 
1 

V COS(PQTO) 



sin(pir L ) \ 

sin (pqt l ) 
1 

cos (pqt l ) J 



(30) 



For arbitrary N , L — 2Q, and detT ^ we can write down the solutions for 
"23b R as 



fc = 0,l,...,L. (31) 



2Q 
k=0 



r„ - [pnl + Im [B n [I - \\b n \\ z b n ®b n 
r„ = Re (s„ (i - \\b n \\- 2 b n ® 6„)) , 



n = 1, Q 
n = 0,1, 



y j j_ , . . . , 



(32) 
(33) 



Since the matrices r„ contain (iV 2 — N^j free parameters, to obtain a unique 
solution one needs to increase the value of Q by factor of (iV 2 — N) to constrain 
the additional degrees of freedom. 
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Some remarks about the meaning of the periodic extensions should be made. 
Since the initial conditions were taken as a linear combination of exponential 
solutions, the existence of the extensions is equivalent to the existence of 2Q + 1 
pair-wise conjugate zeroes of the characteristic function that are located on the 
imaginary axis. As mentioned in the previous section, there can be only a finite 
number of zeroes of the characteristic function in any finite width vertical strip 
of the complex plane. Hence, there could be only a finite number of zeros of 
the characteristic function lying on the imaginary axis. As a result, Q must be 
finite and only finite-dimensional spaces of initial conditions can be extended 
periodically. A practical consequence of this is that one can store a trajectory as 
a periodic extension only approximately, with the error of approximation given 
by the error induced by truncating the Fourier series. However, as follows from 
the theory of Fourier expansions, the error can be made arbitrarily small at 
least in L 2 norm by increasing N. 



4.1.1 Non-linear DDE 

Let us now consider the existence of the periodic extensions in the presence of 
nonlinearities. Take, for example, the extended Hopfield model with dynamics 



L 

X = 

k=0 



J2^k-cr{x(t-T k ))+y{t) (34) 



Proceeding as before we expand the periodic the initial condition in the Fourier 
series for t S [0, Tl] as in ( |2l| ) and obtain that for <p(t) to be a solution of 
homogeneous DDE it needs to satisfy 

L 

4>{t) = Y,A k -a{^{t-T k )) (35) 

fc=0 

Substitution of ( plj ) into ( |35| ) and expansion of the nonlinear term in Fourier 
series yields a system of equations on the unknown matrices A k 



0, n = -Qi, ...,Qi (36) 



ip-nK - ^2 A k ■ ex P (-ipnTk) ■ A„ (n; {b m }) 

V k=0 ) 

where Q\ is not necessarily equal to Q and for each p = 1, ....N the coefficient 
(A„ (r fc ; {b m })) p is defined by 

1 f TL ( Q \ 

A„(T fc ;{6 m }) = — / dte-np(~ip n (t - r fc )) er } exp (ip m (t - T k )) ■ b m \ 

TlJ o \ m =- Q ) 

(37) 

The additional time-independent factor exp (ip n r k ) in the definition of A„ (r k ; {b m }) 
ensures that if a (z) = z then A n (r k \ {b m }) = b n . 
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Now the vector coefficients of the Fourier expansion ( |2l[ ) cannot be factored 
out and the solution for Ak will depend not only on the delays Tk but also on 
6„even for N = 1. Such dependence is desirable in a memory model, since it 
implies that the choice of the parameters Ak is pattern specific. 

When Qi ^ Q and Qi < oo we obtain essentially the same system of equa- 
tions as before. One class of nonlinearities when this occurs is a set that may 
be called polynomial nonlinearities with 

N 2 

a(z) = ^2c n z n , 0<iVi <N 2 <oo (38) 

Ni 

For example, if we assume that all trajectories that we wish to store are bounded, 
in order to introduce a sigmoid nonlinearity that is frequently used in memory 
models, one can take a cubic nonlinearity with 

<t(z) = - (1/3) a 2 z 3 + f3 2 z, a,/3>0 (39) 

and load only such trajectories for which a (z) < (2/3) (/3 3 /a) . Substitution of 
(|9t) into @ results in 



I I i Q W 

ipnK - ^ A k ■ exp (~ip n T k ) ■ f3 2 b n + - 2J b p b q b n - p - q 
V fc =o \ p,q=-Q J J 



0, (40) 



n = -Qi, Qi 

As in the linear case this equation can be still considered as an eigenvalue 
problem for b n , except for a nonlinear operator A = g (A) . On the other hand 
this equation can also be considered as a linear equation on Ak for a given set 
of b n . One novel feature that appears in the nonlinear case is that since we 
have a cubic nonlinearity, in general, up to three different values for each b n are 
possible so that the term (j3 2 b n + | J2pq=-Q b p b q b n - p _ q ^ has the same value 
and hence the solution for Ak is the same. This observation can be used for 
multiple memory storage of a set of M desired memories {b„} ,B = 1, ...,M. 
Consider the following system of equations 



P,9=—Q 



This cubic system of 2Q+1 equations has 3 2< ^ +1 choices of solutions. As a result, 
in general, one can load 3 2( ^ +1 memories corresponding to the same choice 
of matrices Ak- The same argument applies to any polynomial nonlinearity 
defined by (pq) with the corresponding maximum number of memories growing 



as N 2Q+1 . Of course for this storage method to be practical a prescription should 
be given for a way to store a given set of memories for a given nonlinearity, rather 
then choosing a nonlinearity to fit the needed set of memories as we described 



12 



above. In addition a thorough investigation of the basins of attraction of the 
multiple memories needs to be carried out. These and other questions shall be 
addressed in a future publication. 

A problem can arise if Q\ is sufficiently large, that is the system can become 
overdetermined if Q\ is large enough. If Q\ is finite the problem can be cured 
by choosing appropriately large L. The real problem appears when nonlinearity 
is such that Q\ = oo. Then, one would expect that the system ( p5| ) generically 
has no solutions. When Q\ = oo our approach can be nevertheless applied 
approximately in the following sense. If A„ (t^; {&,„}) — > sufficiently fast as 
n — > oo, we can truncate the system at a finite Qmax and the proceed as before 
with computing the Ak- Of course since we zeroed the coefficients A„ (r^; {b m }) 
for n > Q m ax hi effect we substituted the original nonlinearity a (z) with another 
one a (z) and the periodic extensions obtained for a (z) will not be periodic 
solutions for the original nonlinearity. However if the original nonlinear DDE 
is stable with regard to small perturbations of the solutions, the aperiodicity 
will not grow in time and will remain small. Therefore we still be able to store 
time- variable patterns, although with an additional error. 

The criteria for the stability of the periodic extension follow from the sta- 
bility analysis of the associated DDE. For the nonlinear case the stability of 
the nonlinear DDE is equivalent to the stability of the linearized version of the 
original DDE. The basic result about the stability of a linear DDE is that any 
solution of a DDE is stable if and only if the largest real part of the roots of 
the characteristic equation is non-positive. The criteria for the roots of the 
characteristic equation to lie in the complex left half-plane have been studied 
extensively. A summary of the results can be found in . 

5 Parameter Adaptation Dynamics 

If the construction of periodic extensions is realized in the physiology of the 
brain then there must exist an algorithm that adjusts weights and/or delays to 
arrive at the needed values so that a population of neurons learns the appropriate 
weights after many repetitions of exposure to the external input. In this section 
we present examples of such learning algorithms. Of course we do not claim 
that either of the algorithm is actually implemented in the brain. Neither is 
learning a central point of this paper. Nevertheless it is instructive to consider 
some possibilities to estimate the difficulty of the problem. 

We shall consider that the delays are known and fixed and only the weight 
matrices Ak need to be determined. The delays can also be considered as 
variables instead or together with Ak . Investigation of this possibility we leave 
for elsewhere. 

A learning algorithm is not difficult to construct by introducing a Liapunov 
dynamics in the space of weights so that at the fixed point of the evolution we 
obtain the needed periodic extension. The weight dynamics can be consequently 
constructed by defining an error functional that measures the distance from the 
desired solution to the initial condition and using it as a Liapunov function. 
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Consider, for example, "slow" learning done during the multiple presentations - 
"epochs" - of the same initial condition. The "m"th update of A k would then be 
done once per "epoch": one update for each time interval [m tl, {m + 1) tl] ■ 
As an example we can take the "m"th epoch error functional as 

Aw,+1)tl 

V m [A] = / dt \x(t)-cj>o{t)\ 2 , (41) 

J m tl 

where x (t) is the solution of the DDE with a given initial condition 4>o (t) ■ " Fast" 
adaptation within the single epoch can also be defined in a similar way, although 
from the point of view of biology it might be more suitable for adaptation of 
delays. The error functional ( |4l| ) is phenomenologically appealing since both 
x (t) and <j> (t) can be physically " measured" by the encoding population of 
neurons. The discrete evolution in the space of weights with fixed delays can 
then be written as 

{Ak) m+l = (A k ) m+1 - e^-V m [A] , (42) 

where e is a small parameter. If the fundamental solution can be computed 
then x (t) in ( pil"| ) is given by 

x(t)=X(t- tl) 0(t l ) + J2 d T X{t-T- r fe ) A k (t) , 

k=1 Jr L -T k 

where the dependence of the fundamental solution X(t) on the weights is given 
by (0J^). Computation of the appropriate derivatives of V[A] involves computing 
the A— derivatives of the fundamental solution. These can be written out as 



dXpq 

diA^ 



/ dX exp A (t - T k ) (H- 2 (A)) 6 qj , (43) 

Z7Tl J{c) ' 

which enables us to compute gjp^'m [A]. 

A more direct approach to learning the weights is to use the characteristic 
equation itself. We know from the discussion in the preceding sections that 
matrices A k must be chosen in such a way that the characteristic equation has 
purely imaginary roots located at A„ = -^n, n = —Q,..,,Q. Therefore an 
alternative Liapunov function for the weight dynamics can be chosen as 



1 Q 

Vh[A} = - l d et#(A„)| 2 (44) 



-Q 

L 

H{\) = \I-Y,A k exp(-\T k ) 

k=0 



The weight dynamics induced by this Liapunov function is given by 

A k = -eRe\ V exp (-A„r fc ) det H (A„) H~ T (A„) | (45) 
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where bar denotes complex conjugation. This system of equations conceptually 
simpler than the one obtained from ( f4l| ) , since it does not involve integration 
over time. 

Yet another approach is to use the explicit form of the solutions for A k and 
define the norm in the space of matrices by 

1 L 

V Y [A] = -J2tr(A k -Y k )(A k -Y k f (46) 

fc=0 
2Q 

^* = T kn r n 

fc=0 



where T, T arc defined by (2S 2^, 32| [33|) . Differentiation with respect to unknown 
A k results in 

A k = -e((A k -Y k )) 

The three algorithms we presented all apply only to the linear case. For 
applications in biology the DDE are typically nonlinear. In this case one can use 
a stochastic annealing approach to the error functional V m [A] defined by (^) . 
Another alternative is to restrict learning to the linear regime of the nonlinear 
DDE. This is possible because the typical nonlinearities in memory models, e.g. 
the Hopfield model, do have a region of definition of the nonlinearity where 
a (z) w z. 



6 Summary 

In this paper we presented a method for approximate storage of the trajectories 
of nonlinear dynamical systems as periodic or almost periodic solutions of the 
linear and nonlinear delay-differential equations. 

Although the main motivation for this work was to develop a model for time- 
variable pattern storage by the brain, the results might have wider applicability. 
Indeed, they indicate that linear delay differential equations are, in a sense, 
universal models for nonlinear dynamical systems. Choosing appropriately large 
N one in principle can store arbitrary (but finite) number of periodic orbits of 
a nonlinear dynamical system. 

In addition to providing a method for storage of time-dependent patterns our 
model has other features that are attractive for a biological interpretation. For 
example, the method requires global synchronization of populations of neurons, 
since the initial conditions have to be represented via Fourier series with the 
same fundamental frequency ^ . The maximum delay tl can be thought of the 
attention span of the population: periodic trajectories with period larger then 
tl cannot be stored. The model does not have unknown parameters. Once the 
delays are known the weights can be determined uniquely for a given trajectory 
of the environment. 
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There are a number of open questions about our method. Although the 
situation with the linear DDE seems to be clear, a more detailed investigation 
of the effect of the presence of nonlinearity should be carried out. One beneficial 
effect of nonlinearity could be multistability for a chosen set of the coefficients 
Ak,Tk'. the existence of different, initial condition dependent periodic solutions 
similar to the existence of multiple fixed points in the Hopfield model. This 
would make multiple trajectory storage more efficient: a linear DDE can store 
additional orbits only by increasing its dimension linearly with the number of 
the orbits. 

We have shown that nonlinearity brings along one interesting feature: mul- 
tistability. More than one periodic solution can be stored for a given set of 
matrices Ak . When nonlinearity is polynomial, one can estimate that maximum 
possible number of memories grows exponentially in the number of Fourier co- 
efficients. Additional indication that the multistability should exist comes from 
considering a continuum limit of ( |35| ) , a limit that is reasonable to take when 
there are large number of delays that depend smoothly on their index. This is 
the situation when neurons in the brain have large dendritic trees. With obvious 
definitions the continuum limit can be formally written as 



for A a parameter and for some monotonic delay function r (s) . Here we nor- 
malized the matrix- valued function A (s) so that tr (A(s) T A(s)) = 1. With 
the use of the Fourier transform this equation can be related to the Volterra 
and Hammerstein classes of nonlinear integral equations, which exhibit bifur- 
cations in the parameter A. The theory of nonlinear integral equations has 
been exhaustively studied and a number of criteria are available for determin- 
ing the bifurcation points flio| |. These depend mainly on the growth properties 
of the nonlinearity a (z). A more detailed analysis of ( [f7| ) shall be considered 
elsewhere. 
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